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Abstract 

We discuss precision Monte Carlo (MC) calculations for solving the QCD evolu- 
tion equations up to the next-to-leading-order (NLO) level. They employ forward 
Markovian Monte Carlo algorithms, which provide rigorous solutions of the above 
equations. These algorithms are implemented in the form of the Monte Carlo pro- 
gram EvolFMC. This program has been cross-checked with independent, non-MC, 
programs (QCDNuml6 and APCheb33) and the numerical agreement at the level of 
0.1% has been found. 
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1 Introduction 



Evolution equations of the quark and gluon distributions in a hadron, known as the 
DGLAP equations, derived in QED and QCD using the renormahzation group or dia- 
grammatic techniques [1] can be interpreted probabihstically as a Markovian process, see 
e.g. Ref. [2]. Such a process can be modeled using Monte Carlo methods. The correspond- 
ing MC algorithm provides, in principle, an exact solution of the evolution equations for 
parton distribution functions (PDFs). In practice, the main limitation of such a solution 
is the size of a generated MC sample, i.e. corresponding statistical errors of numerical 
results. This is probably the main reason why this possibility has not been exploited 
until recently. Instead, alternative numerical methods and programs solving the QCD 
evolution equations much faster than the Markovian MC have been used, see e.g. [3-5]. 

The feasibility of solving efficiently the DGLAP equations [1] at the leading-order (LO) 
approximation with the Markovian MC was demonstrated for the first time in Ref. [6] . The 
main conclusion of the above work was that the currently available computer CPU power 
allows to solve efficiently and precisely (at the per-mill level) the QCD evolution equations 
with the use of the Markovian MC algorithm. Of course, this method will always be slower 
in CPU time than non-MC techniques. However, it has several advantages, such as: no 
biases and/or numerical instabilities related to finite grids of points, use of quadratures, 
decomposition into finite series of polynomials, accumulation of rounding errors, etc. It 
is also more fiexible in treatment of the PDFs (e.g. no need to split them into singlet and 
non-singlet components) and easier to extend to higher orders, new contributions, etc. 
The above Markovian algorithm can form a basis of a final-state radiation (FSR) parton 
shower MC program, which not only solves numerically the evolution equations but also 
generates events in terms of parton flavours and four-momenta. Moreover, this algorithm 
is a starting point and a testing tool for various kinds of constrained MC algorithms being 
developed for the initial-state radiation (ISR), see e.g. Refs. [7-10]. 

Here we briefly discuss the Markovian MC solution of the DGLAP evolution equations 
up to the next-to- leading order in the perturbative QCD; more details can be found in 
Ref. [11]. The paper is organized as follows. In Section 2 we present a general structure of 
the DGLAP equations and discuss their basic features up to the next-to-next-to-leading 
order (NNLO). In Section 3 we briefly present the Markovian MC algorithm for parton- 
momentum distributions. Numerical results from EvolFMC at the NLO are presented in 
Section 4. They are compared with the results of non-MC program APCheb33. Com- 
parisons with another non-MC program, QCDNuinl6, are also briefly discussed. Finally, 
Section 5 contains the summary and outlook. 
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2 QCD evolution equations 

The general form of the DGLAP evolution equations reads 
d 



din jj/ 
d 

din jj? 
d 

din jj? 



j 
j 

G ^ J2 {^Gg, ^ qj + Pan, ® Qj) + Pgg^G, 



where {qi, . . . , g^^,, q^, . . . , G}{ii, x) - quark, antiquark and gluon distributions; x - 

Bjorken variable; - hard scale, (e.g. = y^Q^ in DIS). 

The integral convolution denoted by (8) involves only longitudinal momentum fractions: 



1 1 



{P ®q){li,x) ^ j dy J dzS{x - zy) P{as,z)q{ii,y) 

\ ° (2) 

/dz / x\ 

— P{as, z) q (^H, -j . 

X 

The splitting functions P{as, z) depend on /i through the coupling constant = CKgi/J')' 
Pia.,z) = I^P^^\z) + (I)'P(^H^) + il^yP'^Hz) + ... . (3) 



— v — 

NLO NNLO 



From the charge conjugation and the SU {uf) symmetry the splitting functions P have 
the following general structure 
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This leads to the basic form of the DGLAP evolution equations 

j 3 

+ Pfg^G, 

^ Qi =Pgg ®Qi + ®Qi + Pm ® + Pm ® Y.^j 



^ 3 3 \P) 

+ Pfg^G, 

3 

+ Pgg®G. 

Within a given approximation some splitting functions may vanish or be equal, e.g. at 
the LO: P;/) = P^f = P,?°) = 0, and at NLO: P^}'^ = P^'^ . 

2.1 Singlet case 

The singlet PDF is defined as 

^ ^[qj{n,x) + qj{n,x)] . (6) 

Introducing the notation 

pff-p: + n,p^ , p7-p^f + > (7) 

we obtain the following evolution equations for the quark-singlet and gluon distributions 

d 



„^ -E = P^^® S + (2n/Pre)®G', 
d 

— — - G ^Pgf®^ + Pgg®G . 
am /j,^ 

The above sphtting functions obey the general relations 
1 

J dz{zPFF{l^:Z) + zPgf{I^,z)} 


1 

= J dz{2nfzPFG{fJ',z) + zPGG{fJ',z)} = 0. 



(8) 



(9) 
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This leads to the momentum sum rule 

1 



J dx {xT,{fj,, x) + xG{fj,,x)} — const, (10) 



where const = 1 in the parton model. 



2.2 Non- singlet case 

The basic non-singlet PDF reads 



= Y] -qj{l^,x)] , (11) 



and its evolution equations is given by 

d 



9 In ij? 

where the new splitting function 



V^P}^s®V, (12) 



PL = P'^ + nfP'_, P^' = P^-P^. (13) 

The set of the splitting functions (the QCD kernels) usually represented in the litera- 
ture reads 

{Pi, Pi, Pfg, Pgf, Pgg}- (14) 

Pf = at the LO, = at the LO and at the NLO, others 7^ at any order. Having 
the above splitting function one can write and solve the evolution equations in any of the 
presented forms. In our Monte Carlo approach we work directly in the flavour space. The 
general parton-parton transition matrix for a gluon and three quark flavours (d, u, s) as 
well as its LO and NLO contributions are given explicilty in Ref. [11]. 

2.3 Behaviour at z ^ 1 

The splitting functions {P± , P^, Pgg} have the following form 

Pias, z) = .f^"^'] + B{as) S{1 - z) + P{a,, z) . (15) 
(1 - ^)+ 

The functions A{as), B{as) and P(a;s, z) are calculated in powers of ctg, e.g. 

fe=0 
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where at the NLO and the NNLO the coefficients D^^\z) are logarithmically divergent: 

= Dk ln(l - + 0{l) . (17) 

Similarly, the splitting functions {PfGi Pgf} contain logarithmically divergent terms: 

r 0{a,) at LO (k = 0) 

P{as,z) ^ I o\al\ii^{l- z)) atNLO(k=l) (18) 
\ 0{alhi\l - z)) at NNLO (k - 2). 

This can lead to big positive or negative weights in Monte Carlo computations. 
2.4 Behaviour at z — ^ 

The splitting functions {-P^ , -P^} are logarithmically divergent at 2; = starting from the 
NLO 

fe=0 i=l ^ 

The remaining sphtting functions {P+, PfG: Pgf, Pgg} have the following behaviour: 

P{as, z) = Ei{as) — + E2{as) - + C(ln''=^) , (20) 

z z 

The logarithmic term is present starting from the NLO {k = 1) approximation: 

E^{a,) = alE^^^+alEf^ + ... , (21) 
while the 1/ z term is present from the LO {k = 0) approximation 

E^{as) = as Ef^ + E^^ + a] Ef\.. . (22) 

3 Markovian MC for parton-momentum distributions 

In Ref. [11] we have described a Markovian MC algorithm for parton distributions and 
we have implemented it in the MC program. However, the factor in the brems- 
strahlung kernels causes a significant loss of MC efficiency! We can get rid of this annoying 
phenomenon by switching to the xD{x) which evolve with the kernels zP{z). The reason 
for improvement is that the kernels zP{z) fulfill the momentum sum rule. 
The evolution equations for xD{x) read 

1 

/dz OC f 3j \ 

— z'S>Kjit,z) -Dj(t,-) . (23) 
z z \ z / 



5 



The kernels z) = 2Pxj{as{t), z) are split into virtual and real contributions: 

z) = t{t)) Skj 6{1 -z) + 7%j{t, z), 

'?%j{t, z) = ^kAI, z) e{\-z- e{t)) e{z - e'), 

where e is an infra-red (IR) cut-off. 

The iterative solution obtained from the above formulae reads 

xDK{t,x) =e-*^(*'*"\TL)j^(to,x) 

1 t 1 

+^ n I ~ ^^-i) / 

n=l f, Ko,...,Kn-i i=l 



(24) 



to 



i=i - 

n 

X Xo-Di^o (^0, 2^o) 5(x - Xo Zi) , 



(25) 



i=l 

where K = Kn- 

The running as{t) can be absorbed into the evolution variable by the transformation 

t 

dt _ as{tA) 



t > T 



With the choice of af'\t) in the definition of r and tA = to we get the iterative solution 
xDk{t,x) =e-*^(^'""°)xL';^(To,x) 

OO j, n 

+E / 5^ n 



1 /" , , , dt CKcl 

Qis(^A) J or a, 

tA 



(26) 



n=l 



i<'0,...,i<'ri-l «=1 

n r 



dTiQiji -Ti-i) I dzi 
( 

i=l - 

n 

X xo-Di<ro('^o,2;o)(^(2^-a:o JJ^i)' 



(27) 



where 



a^'^\to) 



«r(t.) 



(28) 



In order to generate the above distribution with the MC methods we simplify the 
QCD kernels 



(0)/. N 

VfAr. z) - Pf^(ro, z)^e{l-z- e)^^zPi%) , 



zPi%) 



(0) 



.(0), 



(29) 



The approximate kernels do not depend on r! The compensating weight is 



=1 '^KiKi. 

The probabihty of the forward Markovian leap is now 

oo 1 

J dn J dzi ^ uj{Ti, Xi, Ki\Ti^i, Xi^i, Ki_i) = 1 , 

Ti-l 

The real-emission form factor is defined as follows 

n 1 



Ti-l 



n 



= {Ti - Ti-l) ^ TTjK = in - Ti-l) Rk ■ 



(30) 



(31) 



(32) 



On the other hand, the exact virtual (Sudakov) form factor is 



T 

^k{t,to) = j 



TO 



af\t') 



At the LO, for the one-loop af^ and e(r) = e = const, it becomes simply 



AKK{T')\n^-BKK{T') 



^k{t,Tq) = (r-To) 



«f (to) 



TT 



In - - R 



?(o) 



(33) 



(34) 



At the NLO it is much more comphcated, nevertheless it can also be integrated analyti- 
cally, see Ref. [11]. 

To complete the Markovianization, the integral over the "spill-over" variable r^+i is 
added with the help of the identity 



g-*K„(-r,T„) _ gAK„ (T,r„) 



X 



oo 1 

Kn), 

_ n KnJ-1 



(35) 



where Zn+i = x„+i/x„, and 

^Kin, Tj_i) = fxin, Ti_i) - ^Kin, Ti_i) 

The advantage this method is that at the LO for e = e we obtain 

Ak^O, (37) 

due to the fact that the kernels obey the momentum sum rule. This is also valid at the 
NLO in the MS scheme. In the actual MC calculations, Ak can be non-zero due to 
simplifications in the QCD kernels at the low MC generation level. 

The final formula for this MC scenario with the importance sampling for the running 

as reads 



xDKiT,x) 



ri>T ^1 



1 



(38) 



/I /I I ii n 

+ ^ j dXQ j dTn+ldZn+l X] 5Z 11 / ^'^i^^i 
r„+i>T ^"+1 Ko,...,K„-i j=V.<^ 

n 

X (jj{Tn+l,Xn+l, Kn+l\Tn, Xn, K^) JJ^ w(Ti, Xi, Ki\Ti-i, Xi-i, Ki_i) 

1=1 

n 

X S{x - XoY\_ Zi) XoDko{to, Xq) Wp Wa ■ 
i=l 

where Zi — Xi/xi-i, K = Kn and 

n 



i=l 



For explicit expressions of all ingredients of the above formulae and for more details see 
Ref. [11]. 



4 Numerical tests 

We have implemented the above Markovian MC algorithm up to NLO in the MC program 
EvolFMC. Then we have performed comparisons of the MC solution of the DGLAP with 
the solutions provided by the non-MC programs QCDnuml6 [3] and APCheb33 [4]. Wc have 
evolved the singlet PDF for gluons and three doublets of massless quarks from Qq = 1 GeV 
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to Q = 10, 100, 1000 GeV. In our test we have used the following parameterization of the 
starting parton distributions in the proton at Qo = 1 GeV: 



xDg{x) = 1.9083594473 ■ x~°-^{l - xf-°, 

xDq{x) = 0.5 • xDsei,{x) + xD2u{x), 

xDq{x) = 0.5 • xDsc!,{x) + xDd{x), 
xDsea{x) = 0.6733449216 ■ x^^-^{l - xf'^, 
xD2u{x) = 2.1875000000 ■ x°-^(l - 

xDdix) = 1.2304687500 ■ x°-^(l - 
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Figure 1: The upper plot shows the quark distribution xDg{x,Qi) evolved from Qq = 1 GeV 
(black) to Qi = 10 (red), 100 (green) and 1000 (blue) GeV, obtained in the NLO approximation 
from EvolFMC (solid lines) and APCheb33 (dashed lines), while the lower plot shows their ratio. 

In Ref. [11] we have presented the results of the comparisons between EvolFMC and 
QCDnuml6 for the gluon and quark-singlet distributions. The agreement at the level 
of ~ 0.1% has been found for both the LO and NLO evolution equations. Here, in 
Figs. [Hand [2] we show the results of the comparisons between EvolFMC and APCheb33 for 
the NLO evolution. APCheb33 solves the evolution equations with the use of Chebyshev 
polynomials [4]. As one can see, the gluon and quark-singlet distributions from the two 
programs agree within ~ 0.1% (the similar agreement has been found also at the LO). 
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Figure 2: The upper plot shows the gluon distribution xDG{x,Qi) evolved from Qq = 1 GeV 
(black) to Qi = 10 (red), 100 (green) and 1000 (blue) GeV, obtained in the NLO approximation 
from EvolFMC (solid lines) and APCheb33 (dashed lines), while the lower plot shows their ratio. 

5 Summary and outlook 

We have constructed the Markovian Monte Carlo algorithm for solving the QCD DGLAP 
evolution equations at the NLO. We have implemented this algorithm in the MC program 
EvolFMC (in C++). We have cross-checked EvolFMC with the non-MC programs QCDnuml6 
and APCheb33, and found the agreement at the per-mill level. MC computation for the 
NLO evolution is ~ 5 times slower than for the LO evolution. Singular behaviour of the 
NLO PpG and Pgf splitting functions at z ^ 1 leads to large positive weights for the 
F G transitions and to negative weights for the G ^ F transitions in the region of 
z > 0.95. This shows the need for additional resummation in this region. So far only 
massless quarks have been considered, however, adding heavy quarks can be accomplished 
rather easily. Also extension to the NNLO seems to be straightforward. This program can 
be used as a testing tool for constrained MC algorithms for the ISR, see e.g. Refs. [7-10]. 
Last but not the least, this algorithm can form a basis for the FSR parton shower MC 
event generator. 
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